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Summary 


Research  on  constructing  adaptive  grid  generation  schemes  has  been 
supported  for  the  past  42  months  under  AFOSR-85-0195 .  This  research  has 
primarily  centered  on  identifying  ways  of  using  the  familiar  Poisson  grid 
generators  as  adaptive  grid  schemes.  The  control  of  mesh  point  spacing 
has  been  achieved  by  deriving  the  relationship  between  the  arc  length  or 
cell  volumes  and  the  grid  control  functions  of  the  Poisson  grid  generators 
Applications  of  these  schemes  have  been  made  to  a  variety  of  flow  problems 
These  concepts  have  been  applied  to  the  generation  of  both  structured  and 
unstructured  adaptive  grids.  Activities  supported  under  this  grant  are 
summarized  and  details  of  the  status  of  current  research  are  presented. 


Introduction 

The  original  grant  AFOSR-85-0195  was  targeted  toward  applying  the  con¬ 
cept  of  equidistribution  of  a  weight  function  over  arcs  or  cell  volumes. 

Early  in  this  program,  it  became  apparent  that  application  of  strict  equi¬ 
distribution  to  produce  adaptive  grids  was  very  difficult.  Lack  of  smooth¬ 
ness  as  well  as  skewness  created  problems  when  the  grid  scheme  was  coupled 
with  a  numerical  method  for  solving  fluid  flow  problems.  These  difficulties 
persist  in  two  and  three  dimensions  and  relegate  the  application  of  simple 
arc  equidistribution  to  one  dimension.  Applications  in  more  than  one  dimen¬ 
sion  produce  generally  unacceptable  results. 

Research  under  this  program  was  redirected  in  order  to  avoid  the  short¬ 
comings  of  arc  equidistribution  laws.  The  implementation  of  strict  arc 
equidistribution  as  the  control  law  in  adaptive  grid  generation  was  abandoned 
and  was  used  as  a  guide  in  evaluating  the  potential  of  other  approaches. 

Using  this  philosophy,  the  standard  Thompson  grid  generation  scheme  was  writ¬ 
ten  in  a  form  that  can  be  identified  as  a  nonlinear  equidistribution  scheme. 
The  relationship  between  the  grid  control  functions  of  the  Thompson  scheme 
and  the  weight  functions  of  equidistribution  was  analytically  established 
providing  a  mechanism  where  the  existing  Thompson  method  can  be  directly 
used  to  generate  adaptive  grids.  This  development  has  been  adequately  re¬ 
ported  in  the  literature  and  will  not  be  repeated  here. 

The  use  of  Thompson’s  method  as  a  nonlinear  equidistribution  scheme 
permits  control  of  differential  arc  lengths  in  each  computational  coordinate 
direction.  This  requires  that  one  control  or  weight  function  be  specified 
in  each  computational  coordinate  direction.  It  is  advantageous  to  compute 
the  solution  for  an  adaptive  grid  where  area  in  two  dimensions  or  volume 


in  three  dimensions  is  controlled.  This  requires  that  only  one  weight 
function  be  specified  as  opposed  to  two  or  three  when  arc  elements  are 
manipulated.  A  method  where  the  grid  control  functions  were  directly  re¬ 
lated  to  cell  volume  was  developed.  This  method  is  based  on  the  original 
Winslow  diffusion  idea,  one  weight  function  may  be  specified  and  an  adaptive 
grid  controlling  cell  volume  can  be  created.  This  work  has  also  been  re¬ 
ported  in  publications  and  lectures  and  details  are  omitted  here. 

The  identification  of  the  grid  control  functions  with  weight  functions 
of  arc  equidistribution  and  the  diffusion  with  cell  area  in  the  usual 
Poisson  grid  generation  equations  are,  in  the  author's  opinion,  the  most 
significant  contribution  to  grid  generation  technology  for  a  number  of  years. 
For  a  long  time,  grid  point  distributions  were  determined  by  using  certain 
types  of  functions  for  grid  control  without  any  understanding  of  why  arc 
lengths  or  cell  volumes  were  influenced.  The  control  of  arcs  and  volumes 
through  the  grid  control  functions  was  developed  analytically  and  verified 
by  numerical  experiment  under  this  grant.  This  is  a  major  breakthrough  in 
constructing  static  grids  and  understanding  point  location  and  control 
function  selection.  Adaptive  grids  can  now  be  computed  because  the  cause 
and  effect  relationships  have  been  established.  From  a  practical  point  of 
view,  most  government  laboratories  and  industrial  organizations  have  Poisson 
grid  generators  in  house.  With  nearly  trivial  code  changes  and  almost  no 
software  development  cost,  the  existing  code  can  be  modified  and  used  as 
an  adaptive  grid  scheme. 

The  Thompson/Winslow  equations  for  generating  grids  form  an  elliptic 
system.  Some  effort  was  expended  in  evaluating  the  practicality  of  using 
the  parabolic  form  of  these  equations  to  produce  an  adaptive  grid.  The 


solution  of  the  elliptic  system  is  rather  slow  and  imp  .mentation  of  the 
parabolic  approximation  ideally  should  be  much  more  efficient.  Some  suc¬ 
cess  with  this  idea  was  achieved  and  the  ideas  have  been  published. 

The  Poisson  equations  can  be  used  to  generate  arbitrary  mesh  systems 
and  the  control  of  cell  areas/volumes  can  be  achieved  using  the  equidis- 
tribution-grid  control  relationships  developed  under  this  program.  As 
a  result,  effort  has  been  expended  in  constructing  both  structured  and  un¬ 
structured  grids.  This  has  been  successful  and  results  have  been  published 


in  the  literature. 


Research  Status 


4 


f 


The  results  of  research  are  normally  presented  at  conferences  or 
published  in  archive  journals  as  results  are  produced.  Since  this  is  a 
final  report,  a  number  of  topics  are  incomplete  and  papers  have  not  been 
submitted  for  publication.  While  the  status  of  these  efforts  has  been 
reported  earlier,  a  brief  statement  noting  the  current  disposition  of 
the  research  is  necessary. 

1.  Orthogonal  Adaptive  Grids. 

Work  on  structured  two-dimensional  orthogonal  mesh  generators 
has  been  essentially  complete  for  several  months.  It  remains  to 
submit  the  development  for  publication.  The  basic  governing  equa¬ 
tions  for  generating  such  grids  can  be  derived  by  discarding  the 
cross  derivative  in  Thompson's  grid  equations  and  writing 

a  (r  +  4>r  )  +  y  (r  +  iT r  )  =-  0  (1) 

C  nn  n 

where 


r  =  (x,y) 


2^2  ^2 

a  =  X  +  y  =  S 

p  ■'n  n 


2^2  -2 

v  -  *5  +  yj  - 


In  terms  of  the  arc  ratio  for  orthogonal  grids 


r+7^r+-^(r  --^r)=0 

f  C  ^2  nn  f  n 


(2) 


where  S  is  arc  length  and 

S 

f 


The  adaptivity  is  controlled  in  this  system  through  f.  The  role  of  the 
grid  control  functions  is  taken  by  f  and  the  arc  derivatives  serve  as 
the  means  of  exercising  adaptive  control  in  the  grid. 


The  arc  derivatives  are  most  easily  established  by  using  the 
concept  of  equidistribution.  Let  and  be  weight  functions  for 

S  and  S  respectively.  Thus 

t,  T) 

S.w,  -  Cj 

s  w  =  c 
n  2  2 

For  simplicity,  let  Cj^  and  c^  be  unity.  The  arc  ratio,  f,  may  then  be 
written 


f  =  —  O 

"2 

Both  derivatives  of  f  in  the  governing  equation  are  logarithmic 
aiogf  _  1 

ac  w^  "  w^  35 

This  form  for  grid  control  in  boty  5  and  n  leads  to  problems  where 
Wj^  and  w^  are  nearly  equal.  In  that  case,  adaption  literally  vanishes 

because  the  derivative  of  w^  cancels  that  of  w^ . 


The  form  of  the  weight  functions  is  usually 


wi  =  1  4-  A(|f)' 


Success  has  been  achieved  by  letting  A  change  as 


A  =  1  +  b 


3u| I 3u I 

ac'  ‘35' 


where  b  is  a  constant.  The  effect  of  this  form  is  to  prevent  the  re- 

2  2 

duction  of  adaptivity  in  regions  where  u^  ^  Application  of  these 


developments  will  be  included  in  a  paper  soon  to  be  submitted  for 
publication . 

Linear  Elliptic  Adaptive  Grid  Generation 

Several  popular  grid  generation  methods  are  based  on  a  set  of 
linear  elliptic  governing  equations.  However,  these  linear  equa¬ 
tions  are  transformed  from  physical  (x,y)-space  to  computational 
(?>n)-space.  The  transformed  equations  then  become  non-linear,  and 
are  solved  in  computational  space  for  the  (x,y)  coordinates  of  the 
grid  points.  Although  the  original  equations  are  linear  and  un¬ 
coupled,  the  transformed  equations  are  non-linear  and  strongly  coupled 
Therefore,  the  possibility  of  solving  the  original  linear  equations 
in  physical  space  is  attractive  and  could  result  in  the  development 
of  a  new  adaptive  grid  generation  method. 

The  basis  for  this  new  adaptive  grid  generation  method  is  solv¬ 
ing  a  set  of  linear  elliptic  partial  differential  equations  in  physi¬ 
cal  space.  As  such,  it  is  necessary  to  have  an  initial,  base  grid 
on  which  the  solution  is  to  be  found.  This  base  grid  can  be  gener¬ 
ated  by  any  method  as  long  as  it  is  reasonably  smooth  and  non-over- 
lapped.  Once  the  base  grid  is  defined,  grids  can  be  generated  that 
are  smooth,  adapted  to  specified  weight  functions,  and  orthogonal 
at  specified  boundaries  (structured  case).  The  basic  elliptic  equa¬ 
tion  set  is  that  of  the  diffusion  grid  generation  method. 

The  diffusion  method  of  Anderson  has  proven  to  be  a  very  power¬ 
ful  adaptive  grid  generation  method.  The  governing  equations  are 


given  by  the  set 


7  •  DVn  =  0 


where  D(x,y)  is  the  'diffusion'  function.  The  diffusion  function  D 
(x , , '  can  be  shown  to  be  directly  related  to  the  local  cell  area  (i.e 
the  Jacobian)  of  the  grid  that  results  from  solving  the  diffusion 
equation  set.  Therefore,  if  one  wishes  to  have  a  grid  that  adapts 
to  a  positive  weight  function  W(x,y),  then  one  should  chose  the  dif¬ 
fusion  function  such  that 
DW  =  constant 

Thus,  where  the  weight  function  is  large,  the  diffusion  will  be  small 
and  in  turn  the  local  grid  spacing  will  be  small.  There  are  several 
advantages  to  this  method: 

1.  The  method  satisfies  an  extremum  principle  for  any  choice  of 
D(x,y)  ^0.  This  means  that  the  method  is  guaranteed  to  pro¬ 
duce  grids  without  overlap. 

2.  Only  one  grid  control  function  (namely  D)  need  be  chosen.  This 
can  be  contrasted  to  other  elliptic-based  methods  that  require 
the  specification  of  two  grid  control  functions  (e.g.  the 
Thompson  scheme  requires  the  choice  of  two  functions  P  and  Q) . 

3.  The  diffusion  equations  can  be  transformed  into  equations  that 
are  identical  to  the  Thompson  equations  where: 

2  -7D 

7  C  =  P(x,y)  =  •  7^ 

2  -VD 

V  n  =  Q(x,y)  =  •  Vn 

The  normal  practice  when  solving  elliptic  grid  generation  equa- 


(2) 


tions  is  to  interchange  the  dependent  (5,n)  and  independent  (x,y) 
variables  to  get  an  equation  set  of  the  form  (for  the  diffusion  method): 


trivial.  In  many  cases,  it  is  desirable  to  have  the  grid  lines 
be  orthogonal  at  certain  boundaries  of  the  domain.  This 
results  in  a  Neumann  boundary  condition  which  can  be  imposed 


4.  directly  on  the  linear  grid  equations.  This  is 


opposed  to  methods  that  solve  the  equations  in  computational 
space.  In  that  case,  the  grid  points  must  slide  along  the 
boundary  and  therefore  the  boundary  must  be  spline-fit  after 
each  iteration.  This  incurs  additional  expense  and  program 
complexity . 

5.  No  interpolation  is  needed  during  each  iteration  to  update 
the  values  of  the  weight  function  at  the  grid  points.  In  a 
standard  method,  the  weight  function  must  be  updated  after 
each  iteration,  which  is  a  very  costly  procedure. 

The  real  advantage  of  this  method  will  be  experienced  in  the  gener¬ 
ation  of  three-dimensional  adaptive  grids.  Standard  three-dimensional 
adaptive  grid  generation  methods  require  that  the  grid  points  slide 
along  boundaries.  This  requires  a  sophisticated  surface  definition 
method  (such  as  bi-cubic  patches).  On  the  other  hand,  the  new  linear 
adaptive  method  does  not  require  the  movement  of  points  along  any 
boundary. 

The  method  was  tested  on  a  set  of  model  problems  set  in  the  unit 
square.  The  first  example  (Fig.  1  )  shows  the  result  of  adaption  to 
the  function 

T,/  V  1  if  r  <  ^4 

W(x,y)  =  ^  (3 

2  if  r  ^  % 

where 

1  2  ^2 
r  =  (x  -  j)  +  (y  -  y) 

The  grid  adapts  well  to  the  weighting  function  and  is  clearly  uniform 


where  W  is  uniform.  The  second  example  (Fig.  2)  involves  adaption 
to  a  model  of  a  shock  wave  where  the  weighting  function  is  given  by 
W  =  1  +  10|Vul  where 

u  =  tanh(2x  -  y)  (4) 

The  grid  adapts  to  the  model  shock  wave  in  a  reasonably  smooth  manner. 

It  should  be  noted  that  the  grid  was  specified  to  be  orthogonal  to 
all  boundaries. 

The  linear  adaptive  grid  generation  method  has  also  been  tested 
on  realistic  turbomachinery  problems.  In  all  cases  it  has  proven  to 
be  robust  and  efficient  and  should  find  wide  application  in  computational 
fluid  dynamics. 

3.  Generation  of  Unstructured  Triangular  Grids  Using  Elliptic  Partial 
Differential  Equations. 

The  past  few  years  have  seen  a  growing  interest  in  the  use  of  un¬ 
structured  triangular  grids  in  computational  fluid  dynamics.  One  of 
the  major  difficulties  in  the  generation  of  unstructured  grids  is  to 
produce  grids  that  are  smooth,  non-overlapped  and  adaptive  to  a 
specified  weight  function.  Therefore,  a  new  method  has  been  devel¬ 
oped  for  generating  adaptive  unstructured  triangular  grids  based  on 
elliptic  partial  differential  equations  (PDE) .  The  use  of  an  el¬ 
liptic  PDE  ensures  that  the  resulting  grid  is  non-overlapped  and 
smooth,  while  adaption  is  achieved  by  using  Anderson's  diffusion 
method  . 

The  new  adaptive  grid  generation  method  is  based  on  solving  the 
equations  (see  above  discussion  for  more  details  on  the  diffusion  method) 


V  •  DV?  =  0 


V  •  DVn  =  0 


(1) 


These  equations  are  transformed  to  computational  space  to  give  a 


set  of  equations  to  be  solved  for  the  (x,y)-coordinates  of  the  grid: 


ax„  -  2Bx,  +  Yx  =  (ax^  -  Bx  -  (Bx  -  yx  )— 
nn  '  C  n  D  ^  C  n  D 

D.  D 

ay._  -  2By_  +  yy  =  (ay,  -  Bv  ):^  -  (By  -  yy  )— 
sn  nn  'n  D  '^n  D 


where 


2  2 

a  =  X  +  y 

n  n 


B  =  x_x  +  y,y 
f;  n  n 


Y 


For  an  unstructured  grid,  the  computational  space  would  also  be  unstruc¬ 
tured  (as  opposed  to  a  structured  grid  that  produces  a  uniform  grid  in 
computational  space).  Therefore,  instead  of  transforming  to  computa¬ 
tional  space,  the  concept  of  a  master  element  is  employed.  Each 
triangular  element  in  physical  space  is  mapped  to  the  single  triangular 
master  element.  All  computations  are  done  on  the  master  element  and 
then  mapped  back  to  physical  space.  The  governing  equations  are  dis¬ 
cretized  using  a  standard  finite  element  procedure  and  solved  by  simple 
point  relaxation.  The  method  was  tested  on  a  model  problem  set  in  the 
unit  circle.  The  initial  unstructured  grid  is  shown  in  Fig.  3  and  was 
generated  using  a  Delaunay  triangulation  procedure.  Figure  4  shows 

the  result  of  requiring  that  the  grid  adapt  to  the  weight  function 


W(x,y)  = 


10 

J, 

10 


if  r  <  5s 

if  r  >  is 


(2) 


where 


^2  ^2 
r  =  (x  -  y)  +  (y  -  y) 

Note  that  where  the  weight  function  is  large,  the  grid  cells  are  small, 
and  vice  versa.  Also  note  that  the  grid  is  relatively  uniform  where  the 
weight  function  is  uniform. 

This  research  has  shown  that  methods  that  have  been  successfully 
used  to  generate  structured  adaptive  grids  can  be  applied  to  generate 
unstructured  adaptive  grids.  The  new  unstructured  grid  adaption 
method  should  find  wide  applicability  in  the  new  flow  solvers  being 
developed  for  triangular  grids. 
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Figure  1  Linear  adaptive  grid  generation-Example  1 


Figure  2  Linear  adaptive  grid  generatlon-Example  2 


Conclusions 


r’lis  program  has  been  a  successful  one  where  major  constibutions  to 
the  construction  of  adaptive  grid  schemes  have  been  made.  These  contri¬ 
butions  were  made  in  applying  Poisson  grid  generation  schemes  as  solution 
adaptive  mesh  generators.  While  application  of  these  concepts  in  a  dy¬ 
namic  environment  is  possible,  the  major  impact  at  present  is  in  using 
the  adaptive  Poisson  generators  in  static  applications. 

Present  techniques  for  dynamic  adaption  are  based  on  integration  of 
the  equations  governing  the  physical  problem  of  interest  and  then  remesh¬ 
ing  and  refining  the  meshes  a  separate  process.  Better  ways  of  generating 
dynamically  adaptive  grids  must  be  developed.  For  example,  if  the  steady 
grid  law  is  dif ferentialed  to  obtain  the  grid  speeds,  the  grid  and  the 
equations  of  motion  can  be  integrated  in  time  to  any  order  and  no  inter¬ 
polations  of  the  dependent  variables  is  necessary.  Research  is  being 
conducted  in  this  area  and  should  be  continued  in  order  to  better  resolve 
the  issues  in  dynamic  adaption  for  both  structured  and  unstructured  grids. 
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Abstract 

Two  schemes  for  constructing  solution  adaptive  grids  are  reviewed  in  this 
paper.  The  first  method  is  designed  to  control  arc  length  elements  comprising 
the  computational  cell  boundaries  while  the  second  directly  controls  cell  area/ 
volume.  These  techniques  are  both  applied  using  the  popular  Poisson  grid  gen¬ 
eration  equations  by  identifying  the  grid  control  functions  in  terms  of  the 
physical  solution.  A  number  of  examples  are  presented  showing  grids  produced 
with  these  methods.  Adaptive  grid  capability  can  be  added  to  existing  Poisson 
grid  generators  with  very  little  additional  investment. 

1.  Introduction 

The  development  of  techniques  for  constructing  solution  adaptive  grids 
has  received  significant  emphasis  over  the  past  few  years.  The  problem  ad¬ 
dressed  in  this  discussion  of  adaptive  grids  is  to  determine  the  distribution 
of  a  fixed  number  of  grid  points  on  a  physical  domain  in  order  to  improve  the 
quality  of  the  numerical  solution.  In  this  approach,  the  grid  points  are 
moved  during  the  calculations  and  the  point  locations  as  well  as  the  vector 

of  dependent  variables  emerge  as  part  of  the  solution.  Solution  quality  can 

not  be  specifically  defined  since  the  exact  objectives  which  dictate  the  numer¬ 
ical  simulation  are  different  in  each  case.  Usually  some  measure  of  global 

error  is  minimized.  However,  in  many  cases,  simply  improving  resolution  in 

regions  of  rapidly  changing  dependent  variable  is  satisfactory. 

Many  schemes  have  been  put  forv/ard  for  generating  adaptive  grids.  Early 
work  was  based  on  the  concept  of  the  equidistribution  of  weight  function  to 
control  arc  length  and  this  idea  worked  well  in  one  dimension  [1,  2].  In  two 
dimensions  the  idea  of  equidistribution  produced  grids  without  any  smoothness 
measure  and  for  high  adaption  rates,  severe  grid  distortion  was  evident  [3,4,5]. 
The  distortion  issue  was  addressed  in  the  work  of  Nakahashi  and  Diewert  [6] 
where  an  angularity  control  was  incorporated  in  the  formulation  for  control¬ 
ling  rate  of  change  of  inclination  of  cell  sides. 

A  more  elegant  formulation  of  the  adaptive  grid  problem  was  provided  by 
Brackbill  and  Saltzman  [7].  In  this  approach,  a  functional  containing  contri¬ 
butions  from  orthogonality,  smoothness  and  adaption  was  minimized.  The  Euler- 
Lagrange  equations  produced  by  this  formulation  form  an  elliptic  system  for 
the  calculation  of  the  grid  point  locations  in  physical  space.  The  advantage 
to  this  approach  is  that  the  relative  importance  of  each  measure  in  grid  con¬ 
struction  can  be  controlled.  The  concept  has  been  used  in  constructing  grids 
for  a  number  of  interesting  examples  in  fluid  mechanics.  The  disadvantage  is 
that  the  method  is  complex  in  the  sense  that  the  grid  equations  are  complicat¬ 
ed  and  require  large  computational  resources  in  two  dimensions.  The  formula¬ 
tion  in  three  dimensions  appears  to  be  impractical  for  the  same  reasons. 


In  more  recent  work.,  Anderson  [8,9]  and  Anderson  and  Steinbrenner  [10] 
have  devised  a  way  of  using  the  Poisson  grid  generation  equations  [11,12]  to 
construct  adaptive  grids.  The  grid  control  functions  are  related  to  weight 
functions  using  the  concept  of  equidistribution  providing  a  way  of  directly 
constructing  adaptive  grids  with  very  little  change  to  existing  software.  Two 
techniques  have  been  developed.  The  first  is  based  on  control  of  arc  elements 
and  the  second  controls  cell  area/volume .  These  two  formulations  are  reviewed 
in  the  following  sections  and  a  number  of  examples  of  the  application  of  the 
ideas  to  fluid  flow  problems  are  presented. 

2.  The  ARC  Control  Scheme 

Perhaps  the  most  widely  used  grid  generation  scheme  based  on  partial  dif¬ 
ferential  equations  is  that  based  on  the  ideas  originally  proposed  by  Winslow 
[13,14]  and  extended  and  implemented  by  Thompson  et  al.  [11].  The  original 
Winslow  scheme  provides  a  system  of  equations  for  a  two-dimensional  case  of 
the  form 

=  0  ,  V^n  =  0  (1) 

for  the  computational  coordinates  (5,n).  These  equations  assure  that  the  grid 
is  the  smoothest  possible  if  the  smoothness  measure  is 

I  =  //  [(7C)^  +  (Vn)^  dx  dy  (2) 

The  Euler-Lagrange  equations  resulting  from  minimizing  Eq  (2)  are  the  grid  Eqns 
(1).  Adding  source  terms  or  grid  control  functions,  P  and  Q,  to  Eqns  (1)  pro¬ 
vides  a  way  of  controlling  interior  point  placement  other  than  through  the 
boundary  conditions.  The  grid  equations  become 

=  P  ,  V^n  =  Q  (3) 

These  equations  are  usually  written  with  computational  coordinates,  as 

independent  variables  yielding  a  nonlinear  partial  differential  equation  for 
(x,y)  of  the  form 

ar..  -  26r,  +  yr  =  -J^  [Pr.  +  Qr  ]  (4) 

where 


r  =  <x,y)^  , 

2^2 
a  =  X  +  y^ 

(5a) 

n  n 

6  - 

2  ^  2 
,  Y  =  +  y. 

(5b) 

and  the  Jacobian  of  the  mapping  is  defined  by 


J  =  x^y  -  X  y_  (5c) 


The  major  problem  in  using  the  Thompson  scheme  (Eq  (4)),  is  selecting  the 


proper  values  of  the  grid  control  functions.  In  references  8,  9  and  10,  the 
relationship  of  the  grid  control  functions  to  the  weight  functions  of  arc  equi- 
distribution  schemes  is  derived.  In  an  arc  equidistribution  scheme  for  a  con¬ 
stant  n  curve  where  arc  length  in  physical  space  is  denoted  by  S  and  C  varies 
along  the  curve,  we  may  write 


S_W  =  C, 

C.  1 


(6) 


where  W  is  the  weight  function  and  Ci  is  a  constant.  Elements  of  arc  length, 
AS,  are  small  where  W  is  large  and  vice  versa.  If  the  Thompson  Eq  (4)  is 
written  along  constant  C  and  constant  n  curves,  the  resulting  system  of  two 
equations  can  be  identified  as  approximating  arc  equidistribution  laws  (Eq  (6)) 
and  consequently,  the  relationship  between  the  weight  functions  and  the  grid 
control  functions  may  be  identified  as 
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a  3W 

.2  W,  — 
J  1 
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3W„ 
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1  3n 


(7) 


(8) 


In  these  expressions,  Wj ,W2  are  the  weight  functions  controlling  arc  length 
along  n  =  constant  and  ?  =  constant  curves  respectively.  In  three  dimensions 
a  similar  expression  is  obtained  in  the  c  direction.  If  the  values  of  ,  W2 
are  coupled  to  the  solution  of  the  physical  problem,  the  grid  on  the  physical 
domain  can  be  made  to  be  adaptive  and  reflect  the  important  features  of  the 
physical  process. 


3.  The  Area/Volume  Control  Scheme 

In  a  1981  report  [14],  Winslow  suggested  constructing  grids  by  employing 
a  diffusion  law  of  the  form 

V- (DV?)  =  0  ,  7'(DVn)  =  0  (9) 

The  computational  coordinates,  (5,n),  are  controlled  by  the  diffusion,  D,  when 
these  expressions  are  used  as  the  appropriate  grid  generation  set.  In  terms 
of  the  Thompson  grid  control  functions 
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2, 


=  P 


7C-VD 
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_  0  _  _  7n  *70 

'*'  r\ 


(10) 


Since  both  P  and  Q  are  expressed  in  terms  of  a  single  diffusion  term  control¬ 
ling  the  production  of  the  computational  coordinates,  an  obvious  question  is 
to  ask  of  the  D  variable  can  be  in  any  way  related  to  cell  area  in  two  dimen¬ 
sions  or  volume  in  three  dimensions.  If  the  relationship  between  D  and  cell 
volume  can  be  derived,  then  the  grid  control  functions  can  be  used  to  produce 
an  adaptive  grid  by  adjusting  cell  volume  using  a  weight  function  approach. [9] 


If  the  Laplacian  operator  is  applied  to  the  inverse  Jacobian 


(11) 


and  it  is  assumed  that  the  quantity 

n  (C  -  C  )  -  C  (n  -  n  ) 

xy  XX  yy  xy  xx  yy 

is  negligible,  it  may  be  shown  that  the  Jacobian,  J,  is  related  to  the  diffu¬ 
sion  D  by  the  expression 


J  =  CD 


where  C  is  constant.  In  terms  of  the  Thompson  grid  control  functions 

D-  0  D  .  D.  D 

P  =  -  ^  +  A  n  -  A  1 

j2  D  D  ’  ~  _j2  D  ■  ^2  D 


(12) 


(13) 


The  actual  adaption  can  be  accomplished  by  assuming  that 

DW  =  constant  (14) 

where  W  is  a  weight  function  related  to  the  solution.  Again,  where  W  is  large, 
D  is  small  (cell  area)  and  vice  versa. 

4.  Examples  and  Conclusion 

It  is  worthwhile  to  show  some  examples  of  grids  computed  using  the  ideas 
presented  above.  Figure  1  shows  the  grid  produced  using  Eqns  4,  7  and  8  for 
supersonic  flow  through  a  compressor.  The  inlet  Mach  no.  was  2.85  and  the 
grid  was  adapted  using  a  weight  function  dependent  on  the  pressure  gradient. 

The  Euler  equations  were  solved  and  the  adapted  grid  shows  the  presence  of  the 
bow  shock  as  well  as  the  passage  shock  between  the  blades.  Obviously,  this 
example  was  computed  with  a  very  course  grid. 

Figure  2  shows  the  grid  used  to  compute  solutions  of  the  Euler  equations 
for  Mach  10  flow  over  a  six  degree  sharp  cone  with  drag  flaps  [15]  in  a  plane 
normal  to  the  cone  axis.  The  form  of  the  weight  function  was  selected  to 
depend  both  on  first  derivature  and  curvature  of  the  pressure  along  computa¬ 
tional  coordinate  lines.  The  nearly  conical  shock  produced  by  the  vehicle  is 
clearly  visible  in  the  clustering  between  the  body  and  the  outer  boundary  and 
the  strong  shock  produced  by  the  drag  flap  is  clearly  visible.  Since  this  is 
a  solution  based  on  pressure  information  in  an  inviscid  flow,  no  other  flow 
details  would  influence  the  grid  point  locations. 

Figure  3  shows  the  grid  point  distribution  employed  to  solve  the  Euler 
equations  for  a  three  dimensional  wing  (23012) .  The  root  section  is  nearest 
the  viewer  and  the  tip  is  into  the  figure.  Not  all  grid  lines  for  this  tran¬ 
sonic  calculation  (M  =  .8,  a  =  2®)  are  shown.  The  grid  was  adapted  using 
pressure  gradient  information  with  Eqns  (7,8).  The  shock  on  the  upper  surface 
is  clearly  visible  and  clustering  toward  the  wing  tip  shows  althougn  the  per¬ 
spective  does  not  provide  a  very  good  view  of  this. 


As  a  final  example.  Figure  4  shows  the  grid  used  to  compute  a  solution 


of  Che  Navier-Scokes  equations  for  normal  injection  in  a  supersonic  stream. 

The  equations  for  laminar  flow  were  used  and  the  weight  function  was  as¬ 
sumed  to  depend  on  the  pressure  gradient  and  the  velocity  gradient.  The 
shocks  and  shear  layers  for  this  example  are  viable  in  the  grid.  In  the  in¬ 
terests  of  brevity,  the  contour  plots  of  flow  variables  have  been  omitted 
but  the  appropriate  references  can  be  consulted. 

Both  the  arc  control  and  cell  volume  control  methods  work  well  within 
the  constraints  imposed  by  the  Poisson  grid  generation  equations.  The  tech¬ 
niques  are  easily  implemented  in  existing  codes  for  a  minimum  investment  of 
time.  Either  procedure  may  be  used  to  rapidly  gain  an  adaptive  grid  capabilic 
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ABSTRACT 

A  technique  for  generating  solution  adaptive  grids  controlling  cell 
volume  is  presented.  The  grid  control  functions  of  the  popular  Poisson 
grid  generators  are  defined  to  control  cell  size.  An  analytic  derivation 
showing  the  relationship  of  the  grid  control  functions  and  the  cell  vol¬ 
ume  is  presented  and  several  examples  are  given  illustrating  the  idea. 
Only  one  function  must  be  prescribed  to  control  cell  size  when  this  tech¬ 
nique  is  implemented  in  a  Poisson  grid  generator. 

INTRODUCTION 

One  of  the  first  steps  in  applying  a  numerical  method  to  solve  a 
partial  differential  equation  is  to  construct  a  suitable  mesh  or  grid 
for  the  physical  domain.  There  are  many  techniques  for  doing  this  but 
the  most  popular  method  is  that  of  Thompson  et  al  (1).  This  scheme 
employs  the  original  method  proposed  by  Winslow  (2)  but  includes  addi¬ 
tional  control  over  placement  of  interior  grid  points.  This  has  led  to 
the  popular  Poisson  mesh  generators  for  constructing  grids  for  arbitrar¬ 
ily  shaped  regions. 

Recently,  a  significant  interest  has  surfaced  in  generating  adaptive 
grids.  When  an  adaptive  method  is  used,  the  mesh  point  locations  change 
while  the  solution  of  the  physical  problem  is  computed.  Both  the  grid 
point  locations  and  the  original  physical  unknowns  evolve  as  part  of  the 
problem  solution.  The  grid  points  are  moved  during  the  calculation  in 
order  to  improve  the  quality  of  the  result.  This  is  usually  assumed  to 
mean  that  some  measure  of  total  error  in  the  solution  has  been  reduced 
or  physical  resolution  has  been  improved  in  regions  where  large  changes 
in  dependent  variables  occur.  The  issue  of  improved  quality  is  problem 
dependent  and  additional  research  will  be  required  to  adequately  define 
both  grid  and  solution  quality. 

Several  schemes  for  generating  adaptive  grids  have  been  developed 
using  a  variety  of  approaches.  These  include  strict  equidistribution 
schemes  employed  by  Dwyer  et  al  (3),  Gnoffo  (4)  and  Rai  and  Anderson  (5). 
A  modification  of  the  strict  arc  equidistribution  scheme  to  reduce  grid 
skewness  was  used  by  Nakahashi  and  Deiwert  (6)  and  has  been  applied  to 
a  number  of  problems. 


A  number  of  researchers  have  proposed  mesh  schemes  based  on  deriving 
the  grid  generation  equations  by  minimizing  some  grid  function.  Most 
notable  of  these  methods  is  that  of  Brackbill  and  Saltzman  (7)  where  a 
combination  of  contributions  from  smoothness,  orthogonality  and  adaption 
were  minimized.  The  Euler  Lagrange  equations  that  result  provide  the 
elliptic  grid  generating  set  for  calculations  of  the  grid  point  locations. 
As  an  adaptive  mesh  scheme,  this  technique  can  be  applied  and  has  been 
shown  to  work  well  in  some  applications.  The  expense  of  solving  these 
equations  along  with  the  complexity  of  the  system  (especially  in  three 
dimensions)  has  prevented  this  technique  from  being  widely  used. 

More  recently,  Anderson  (8)  and  Anderson  and  Steinbrenner  (9)  have 
used  the  Thompson  scheme  as  an  adaptive  grid  generator.  The  grid  control 
functions  are  related  to  the  solution  of  the  physical  problem  resulting 
in  a  direct  way  of  controlling  point  locations  in  physical  space.  As 
the  physical  solution  evolves,  the  grid  can  be  recomputed  providing  grid 
point  placement  where  desired.  This  represents  a  significant  advance 
because  very  little  software  modification  is  required  to  make  existing 
Poisson  grid  generation  codes  fully  adaptive. 

In  the  following  sections,  a  technique  for  creating  adaptive  grids 
controlling  cell  volume  is  presented.  This  method  is  developed  to  work 
in  conjunction  with  the  Thompson  scheme.  The  relationship  between  cell 
volume  and  the  grid  control  functions  is  based  on  solution  of  a  differen¬ 
tial  equation  relating  cell  volume  and  the  computational  coordinates. 

A  number  of  simple  examples  are  given  illustrating  the  application  of  the 
method  in  two  dimensions. 


A  CELL  VOLUME  ADAPTIVE  GRID  SCHEME 

The  modification  of  the  Thompson  scheme  to  produce  adaptive  grids 
in  Refe  (8)  and  (9)  works  well  and  is  attractive  due  to  its  simplicity. 
However,  the  grid  control  functions  are  modified  to  produce  grids  where 
arc  equidistribution  laws  are  approximately  obeyed.  Consequently,  a 
weight  function  is  required  for  each  dimension.  In  three  dimensions, 
these  weight  functions  are  required  to  provide  the  grid  adaption.  It 
may  be  advantageous  to  prescribe  only  one  control  or  weight  function  to 
adjust  the  cell  volume  rather  than  three  functions  to  independently 
adjust  arc  lengths. 

A  technique  for  producing  a  grid  generation  scheme  controlling  cell 
volume  where  only  one  function  must  be  specified  is  possible.  Winslow 
(2)  suggested  writing  the  grid  generation  equations  in  the  form 


V- (DVC)  =  0 

(la) 

v-(DVn)  =  0 

(lb) 

where  (^,n)  are  the  computational  coordinates  and  D  is  a  diffusion.  The 
diffusion  controls  the  flux  of  computational  coordinate  through  the  sur¬ 
face  of  an  arbitrary  volume.  If  Eq  (la)  is  integrated  throughout  an 
irbitrary  volume  and  the  divergence  theorem  is  applied. 


//DVC-NdA  =  0 
S 


(2) 


From  this  expression,  it  is  clear  that  D  does  indeed  control  the  pro¬ 
duction  of  C  (or  n)  in  the  volume.  This  may  also  be  seen  by  expanding 
Eqs  (la)  and  (lb)  to  obtain  the  more  familiar  form 

2  VD 

^  ?  =  -  VC-  (3a) 

7D 

V“n  =  -  Vn-  —  (3b) 

The  right  sides  of  Eqs  (3a)  and  (3b)  are  recognized  as  the  grid  control 
functions  in  Thompson's  method 

VD 

P  =  -  VC-  ^  (4a) 

VT5 

Q  =  -  Vq-  ^  (4b) 

In  earlier  adaptive  grid  schemes  based  on  Poisson  grid  generators  (8), 

P  and  Q  are  each  explicitly  related  to  the  changes  in  arc  length  in  the 
C  and  n  directions  respectively.  In  the  Winslow  diffusion  formulation, 
the  P  and  Q  control  functions  both  depend  on  a  single  function,  D,  which 
determines  the  point  distribution. 

The  diffusion  term  in  the  Winslow  formulation  can  be  directly  re¬ 
lated  to  the  cell  area  or  volume.  Consider  the  Jacobian  of  the  mapping 
defined  as 

I  =  C  n  -  ?  n  =  (5) 

X  y  y  X  3(x,y) 

2 

for  the  two-dimensional  case.  Since  the  diffusion  is  related  to  V  C  and 
^  h,  a  reasonable  step  is  to  form  V^I. 

V^I  =  v^[c  n  -  C  n  ]  (6) 

X  y  y  X 

Expending,  an  expression  of  the  form 

V^I  =  IV^D  -  VD  •  VI  +  R  (7) 


is  obtained  where  R  is  defined  as 


2  ’  <17)  -  ’  (|f)  • '  <|j>i 


and  D  =  in  D.  Introduce  the  inverse  Jacobian 


J 


(9) 


I 

1 


and  let 


J  =  2.nJ 

This  reduces  Eq  (7)  to 

(J  -  D)  -  •  7(J  -  5)  =  -  JR  '  (10) 

The  term  on  the  right  side,  -JR,  is  assumed  to  be  small  by  comparison 
with  terms  on  the  left.  For  those  cases  where  this  assumption  is  valid 

(J  -  D)  -  7j  •  V(j  -  D)  =  0  (10a) 

and  a  solution  may  readily  be  written  as 

V(j  -  D)  =  0 


or 

J  -  D  =  Constant 

and  D=CJ  (11) 

where  C  is  a  constant.  The  diffusion  term  in  Winslow's  scheme  is  simply 

the  local  cell  area  or  volume  scaled  by  a  constant.  With  this  result, 

the  original  grid  control  functions  may  be  written  in  terms  of  the  diffusion. 


The  original  grid  generation  equations  with  (x,y)  as  independent 
variables  are  now  written 


-  23Xj.^  +  Yx^^  =  [dtx^  “  Sx^l  _5  "  [6x^  -  YX„  ]_n 


nn 


r  ■* 

5  n  p- 


“>'55  ■  -  SyJ  _L  "  [Sy,  -  YyJ_iI 


5n  ^nn 


(12a) 

(12b) 


where 


a 

6 


2 

X  +  y 

n 


x^x  + 

5  n 


2 

n 


The  local  cell  area  may  be  controlled  by  prescribing  the  diffusion,  D, 


where  it  appears  on  the  right  side  of  Eqs  (12a)  and  (12b).  If  an 
adaptive  grid  generation  technique  is  desired,  the  diffusion  or  local 
area  may  be  related  to  the  computed  dependent  variables  of  the  physical 
problem. 

EXAMPLES 

To  illustrate  the  control  of  cell  area  with  this  diffusion  formu¬ 
lation  of  the  Poisson  grid  generation  equations,  a  few  examples  are 
necessary.  Suppose  an  analytic  function  u(x,y)  is  prescribed  on  the’ 
unit  square.  The  grid  point  distribution  can  be  controlled  by  using  an 
area  equidistribution  law.  Define  a  weight  function  W(x,y)  so  that 

DW  =  constant  (13) 

When  the  weight  function  is  large,  the  diffusion  (area)  is  small  and 
vice  versa.  For  the  simple  analytic  examples  presented  here,  it  is  suf¬ 
ficient  to  relate  the  weight  function  to  the  gradient  of  u(x,y). 

W(x,y)  =  1  +  A(7u)^  (14) 

In  this  expression,  A  is  a  constant  although  it  may  be  a  prescribed 
function,  where  the  degree  of  adaption  is  adjusted  by  increasing  or 
decreasing  A. 

The  first  example  is  based  on  adapting  a  31  x  31  point  mesh  to 
the  pillbox  function 

u  =  1.  0<r<.15 

u  =  l.-5r  .15<r<.35 

u  =  0.  .35<r 

The  results  of  computing  a  solution  for  the  grid  point  locations  for 
this  pillbox  function  using  the  weight  functions  given  in  Eq  (14)  are 
shown  in  Fig  1.  Since  the  weight  function  only  increases  in  regions 
where  the  gradient  is  nonzero,  point  clustering  occurs  in  the  region 
where  .15<r<.35.  In  the  regions  where  the  function  is  constant,  the 
grid  is  relatively  uniform.  The  boundary  points  are  allowed  to  move 
along  constant  x  or  y  lines  to  adjust  to  the  adaption  provided  through 
the  grid  control  functions.  On  the  boundary,  it  has  been  assumed  that 
the  second  derivative  of  the  floating  coordinate  with  respect  to  the 
opposite  family  computational  coordinate  is  zero.  For  example,  along 
X  =  0,  y  is  computed  on  a  5  =  constant  surface.  Thus,  is  assumed 
to  vanish  along  this  boundary. 

Fig  2  shows  the  results  of  applying  the  grid  area  adaption  scheme 
to  a  problem  of  a  simple  shock-like  function.  The  prescribed  function 
is  0  to  the  left,  one  to  the  right  and  linearly  transitions  from  0  to  1 
over  15%  of  the  physical  domain  length.  This  example  again  employs  the 
weight  function  of  Eq  (14)  and  illustrates  the  clustering  of  grid  points 


in  regions  of  nonzero  gradients. 

The  grid  control  functions  of  the  Thompson  scheme  can  be  ex¬ 
pressed  in  terms  of  Winslow's  diffusion  and  adaptive  grids  control¬ 
ling  cell  area/volume  can  be  directly  generated  as  evidenced  by  the 
examples  shown.  In  Ref  (8),  arc  length  control  was  employed  to  adjust 
lengths  and  create  a  similar  adaption.  In  cases  where  a  single  com¬ 
putational  coordinate  can  be  made  to  nearly  align  along  an  adaption 
surface,  it  is  best  to  employ  arc  adaption  rather  than  attempt  to  use 
the  volume  method.  However,  in  problems  where  the  form  of  the  solution 
is  unknown,  the  area/volume  technique  has  the  advantage  of  only  requir¬ 
ing  specification  of  one  control  function  through  the  diffusion  term. 

The  linear  relationship  between  the  Jacobian  and  the  diffusion 
given  by  Eq  (11)  is  approximate.  The  term  on  the  right  side  of  Eq  (10) 
has  been  assumed  to  be  small  and  has  been  neglected  in  this  analysis. 

The  meaning  of  this  term  is  difficult  to  establish  and  the  effect  on 
the  diffusion  area  relationship  is  hard  to  assess.  The  examples  pre¬ 
sented  indicate  that  the  influence  is  small.  However,  there  may  be 
cases  where  a  more  careful  evaluation  of  the  impact  of  this  term  is 
necessary.  Any  approximate  analysis  must  be  employed  keeping  the 
restrictions  of  the  method  in  mind. 

SUMMARY  AMD  CONCULSIONS 

The  diffusion  formulation  of  Winslow  has  been  reviewed  as  a  Poisson 
grid  generation  technique.  The  diffusion  was  shown  to  approximately  be 
proportional  to  local  cell  volume  or  area.  The  grid  control  functions 
of  the  standard  Poisson  grid  generators  were  related  to  the  diffusion 
and  several  examples  illustrating  control  of  area  were  presented. 

This  technique  is  another  way  of  using  standard  Poisson  grid 
generators  for  creating  adaptive  grids.  Either  area/volume  or  arc 
length  along  constant  computational  coordinate  surfaces  may  be  control¬ 
led.  The  selection  process  depends  on  the  particular  application  and 
the  skill  of  the  user  in  employing  appropriate  weight  functions  pro¬ 
viding  grid  control.  Like  the  arc  adaption  scheme,  this  volume /area 
control  method  can  be  implemented  with  only  small  changes  in  existing 
Poisson  grid  generators.  The  addition  cf  thrs  option  provides  users 
with  an  additional  choice  in  constructing  useful  grids. 
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ABSTRACT 

Equidistribution  of  a  weight  function  over  a  mesh  is  the  main  concept  employed  in 
the  recent  development  of  adaptive  grid  schemes.  This  idea  is  reviewed  in  this  paper, 
and  several  examples  of  grids  produced  using  equidistribution  are  presented.  The 
shortcomings  of  this  approach  are  identified,  and  its  usefulness  is  evaluated.  Next,  the 
traditional  Poisson  grid  generators  are  written  as  nonlinear  equidistribution  schemes. 
In  this  form,  the  relationship  between  the  numerical  solution  of  a  physical  problem, 
the  source  terms  in  the  grid  generation  equations,  and  the  weight  functions  in 
equidistribution  schemes  are  identified.  Examples  are  given  illustrating  the  ease  of 
providing  a  useful  adaptive  grid  with  this  approach. 


INTRODUCTION 

The  idea  of  constructing  an  adaptive  grid  is  appealing  when  systems  of 
partial  differential  equations  are  solved  using  numerical  techniques.  Concep¬ 
tually,  an  adaptive  grid  can  be  constructed  when  any  type  of  partial 
differential  equation  is  solved.  By  the  very  nature  of  an  adaptive  grid, 
intermediate  results  are  required.  For  hyperbolic  or  parabolic  equations,  a 
marching  direction  exists  and  the  necessary  data  are  available  at  the  end  of 
each  time  or  timelike  step.  For  elliptic  equations,  a  relaxation  scheme  can  be 
used  to  provide  the  evolving  solution.  By  using  the  intermediate  resxolts,  the 
grid  can  be  adjusted  to  position  points  where  a  need  exists  to  refine  the  mesh. 
Normally,  grid-point  motion  is  desired  to  provide  better  resolution  in  regions 
where  rapid  changes  in  the  solution  occur  or  to  provide  some  reduction  in 
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global  error.  These  two  reasons  for  repositioning  grid  points  are  not  mutually 
exclusive,  and  an  improvement  in  one  results  in  an  improvement  in  the  other. 
The  basic  idea  reduces  to  providing  automatic  grid  adjustment  to  resolve 
details  of  the  solution  in  regions  of  varying  length  scale. 

Dwyer  et  al.  [1]  have  computed  solutions  for  several  one-dimensional 
fluid-dynamics  problems  including  adaptive  gridding.  Mesh  points  were 
positioned  in  regions  where  solution  gradients  were  large,  using  the  idea  of 
equidistribution.  Rai  and  Anderson  [2,  3]  employed  an  attraction  scheme  for 
an  adaptive  grid  based  on  equidistribution  of  error  on  the  mesh.  They  solved 
a  number  of  problems  in  both  one  and  two  dimensions  and  illustrated  the 
usefulness  of  a  dynamically  adjusting  grid  in  multidimensional  applications. 
Gnoffo  [4]  also  developed  an  adaptive  grid  scheme  similar  to  that  employed 
in  Reference  [1].  Gnoffo  successfully  applied  his  technique  to  compute 
solutions  of  a  number  of  complicated  fluid-dynamics  problems.  In  one 
dimension,  aH  of  these  methods  work  weU.  Anderson  [5]  unified  many  of 
these  schemes  and  showed  that  the  idea  of  equidistribution  was  common  to 
most  of  them.  A  two-dimensional  adaptive  grid  scheme  based  on  equidistribu¬ 
tion  was  also  presented  in  Reference  [5]  and  was  used  to  demonstrate  that 
lack  of  smoothness  and  control  over  grid  skewing  are  limiting  factors  in  the 
application  of  such  schemes  in  more  than  one  dimension. 

Nakahashi  and  Deiwert  [6]  have  appUed  the  equidistribution  concept  in 
producing  solution-adaptive  grids.  In  addition  to  the  usual  weight  function 
controlling  grid-point  spacing,  grid  skewness  control  was  incorporated  through 
torsion  spring  coefficients.  Results  produced  and  included  in  Reference  [6] 
represent  the  successful  application  of  the  technique  in  both  two  and  three 
dimensions.  However,  the  method  does  not  include  any  smoothness  measure 
and  requires  three  torsion  constants  in  addition  to  the  weight  functions  at 
each  point.  For  these  reasons  alternative  techniques  may  be  desirable. 

Brackbill  and  Saltzman  [7]  have  addressed  the  difficult  multidimensional 
problem  through  a  variational  approach.  A  functional  composed  of  contribu¬ 
tions  from  an  orthogonality  condition,  a  smoothness  measure,  and  adaptivity 
was  minimized,  resulting  in  a  system  of  elliptic  partial  differential  equations 
to  solve  for  the  grid-point  locations.  Mathematically,  this  approach  is  very 
attractive,  and  it  produces  excellent  results  for  physical  problems.  Practically, 
the  method  becomes  nearly  intractable  in  three  dimensions  when  the  grid  is 
calculated  at  each  iteration  level.  _ _ 

It  is  clear  that  a  new  way  of  producing  an  adaptive  grid  which  retains 
some  of  the  properties  incorporated  in  Brackbill  and  Saltzman’s  formulation 
is  needed.  This  scheme  should  require  relatively  modest  amounts  of  computer 
storage  and  CPU  time  for  practical  problems.  Grid  generation  schemes,  such 
as  that  of  Thompson  et  al.  [8],  based  on  solving  a  system  of  Poisson  equations 
are  most  attractive.  The  governing  partial  differential  equations  for  the  grid 
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possess  all  of  the  requisite  properties  of  smoothness  and  relatively  low 
skewness.  The  problem  remains  to  make  this  method  adaptive  by  relating  the 
source  functions  to  the  solution  of  some  physical  problem. 

In  the  following  sections,  the  concept  of  equidistribution  along  arcs  is 
reviewed,  and  results  of  generating  grids  using  strict  equidistribution  for 
several  examples  are  presented.  The  Thompson  scheme  is  then  shown  to 
represent  a  nonlinear  equidistribution  law  when  the  Poisson  grid  equations 
are  written  along  constant  computational  coordinate  surfaces.  Finally,  the 
interpretation  of  Thompson’s  method  as  a  nonlinear  equidistribution  law 
provides  a  direct  relationship  between  the  source  terms  of  the  Poisson 
equations  and  the  solution  of  the  physical  problem  under  consideration.  A 
number  of  simple  examples  are  given  showing  the  utility  of  using  Thompson’s 
method  as  an  adaptive  grid  scheme. 

EQUIDISTRIBUTION 

The  equidistribution  concept  for  generating  a  grid  is  most  easily  explained 
if  a  single  space  dimension  is  considered.  In  this  case,  the  governing  partial 
differential  equation  is  transformed  from  physical  space  (i,  t)  to  a  computa¬ 
tional  domain  (|,  t),  where  the  numerical  solution  is  calculated.  When 
equidistribution  is  used  to  establish  the  mesh,  the  spatial  locations  of  the  grid 
points  are  determined  by  the  equation 

x^w  —  constant,  (l) 

where  w  is  the  weight  function.  With  this  formiJation,  the  weight  hmction  is 
said  to  be  equidistributed  over  the  mesh.  A  new  grid  position  in  the  physical 
plane  is  established  for  each  unit  change  in  This  may  be  seen  by  assuming 
that  A  I,  the  spacing  in  the  computational  domain,  is  set  at  a  value  of  imity. 
With  this  assumption,  (1)  may  be  approximated  by 

Axw  =  c  =  constant.  (2) 

The  physical  mesh  increment  Ax  is  small  when  w  is  large  and  vice  versa. 

The  grid  law  provided  by  (1)  can  be  used  to  generate  a  grid  by  two 
different  methods.  Since  most  researchers  are  familiar  with  second-order 
partial  differential  equations  as  grid  generators,  it  is  appropriate  to  write  (1) 
in  the  form 
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This  expression  can  be  solved  for  the  new  mesh-point  locations  using  stan¬ 
dard  finite-difference  techniques.  Equation  (1)  can  also  be  solved  for  either  x 
or  ^  using  simple  numerical  quadrature.  Let  x  =  at  ^  =  ln,ax>  ^^.d 
integrate  to  obtain 


It  is  worthwhile  noting  that  this  is  the  form  of  the  equidistribution  law 
obtained  when  the  physical  coordinates  are  directly  computed.  The  computa¬ 
tional  coordinate  |  can  be  computed  as  an  alternative  means  of  finding  the 
grid  points.  The  computational  coordinates  are  given  by 


When  the  computational  coordinates  are  computed  with  (5),  the  physical 
coordinates  are  obtained  by  interpolation. 

As  noted  earlier,  the  one-dimensional  equidistribution  law  has  been  used 
successfully  in  computing  solutions  to  a  number  of  physical  problems.  The 
weight  function  w  is  related  to  the  solution  and  typically  is  of  the  form 


lu  =  1  -I-  A 


h 


h 


max 


where  h  represents  a  gradient  or  some  similar  measure  of  the  solution  where 
better  resolution  is  desired,  and  A  is  adjusted  to  alter  the  adaptivity  in  the 
mesh.  The  unity  term  is  inserted  to  prevent  excessive  mesh  stretching  where 
the  value  of  h  vanishes. 

The  equidistribution  law  (1)  has  also  been  written  for  curvilinear  arcs  in 
the  form 


5^10^  =  constant. 


(7) 


where  S  is  arc  length  and  tCi  is  again  a  weight  function.  Equation  (7)  is 
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written  along  a  constant-7j  surface,  and  a  similar  expression  may  be  written 
along  a  constant-^  curve, 


N^W2  =  constant, 

(8) 

where  N  is  arc  length  and  lUg  is  a  second  weight  function. 

As  in  the  one-dimensional  case,  (7)  and  (8)  can  be  solved  for  {x,y)  by 
simple  quadrature  or  by  computing  a  solution  of  the  governing  partial 
differential  ern»ations.  If  the  quadrature  approach  is  selected  Equation 

(4)],  the  changes  in  arc  lengths,  AS  and  AN,  are  computed  and  these  are 
related  to  changes  in  x  and  y  through  the  expression 

Ai  =  ^  AS  +  —  AJV, 

(9) 

Ay  =  ^  AS  +  ^  AN. 

(10) 

New  values  of  (i,  y)  are  determined  by  adding  these  incremental  changes  to 
the  original  values. 

If  (7)  and  (8)  are  differentiated,  the  second-order  expressions  correspond¬ 
ing  to  (3)  are  obtained: 

(11) 

n;, +N-^  =  o. 

W2 

(12) 

In  these  equations 

(13a) 

(13b) 

The  partial  differential  equations  for  (x,  y)  directly  are 

(14) 

XX  +vi/  +fx^  +  !/^)  — ~  =  0. 

W2 

(15) 
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1 1.0,  0  <  X  <  0.3(!/ + 10), 

u  =  (  0-5(7  — 1)  + 0.25 y,  0.5(1/  + 10)  <  x  $  0.5(y  + 14), 

\  0.0,  0.5(t/ + 14)  <  I  ^  24. 


Equations  (14)  and  (15)  form  an  elliptic  system  which  can  be  solved  for  (x,  y) 
when  the  weight  functions  are  prescribed.  Typical  examples  of  grids  pro¬ 
duced  with  this  approach  are  shown  in  Figures  1-4.  In  these  calculations,  the 
weight  function  in  both  the  |  and  tj  directions  is  based  upon  gradients  of 
some  prescribed  function  in  that  direction.  For  example,  if  u  is  a  given 
function  of  (x,  y),  the  weight  fimction  Wi  is  written 


uJi  =  1  -f  A 


'  du 


2 


(16) 


A  similar  expression  for  W2  is  used,  incorporating  the  gradient  in  the  rj 
direction.  The  value  of  u  is  simply  a  prescribed  function  of  (x,  y)  as  noted  in 
the  figure  captions.  With  this  approach,  a  grid  can  be  constructed  with  an 
adaptive  scheme  without  the  complications  of  an  additional  algorithm  to 
compute  the  solution  of  some  other  system  of  partial  differential  equations. 

Several  observations  can  be  made  on  the  basis  of  the  example  problems 
shown.  First,  the  idea  of  equidistribution  can  be  used  to  generate  grids. 
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Fig.  2. 
function. 


Grid  for  circular  notch  problem,  24X24  grid;  equidistribution  with  unity  weight 


Fig.  3.  Grid  for  circular  notch  problem:  Thompson’s  scheme,  P,Q  =  0. 
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Fic.  4.  Equidistribution  grid  for  sine-wave  problem: 


/  0,  0  ^  t/ ^  ll-r4sin(irc/12), 

0.5[!/ -  11  -  4sin(fl’x/12)],  ll-i-4sin(irx/12)  «  t/ «  13  +  4sm(wi/12), 
\  1.0,  13'4-4sin(fl’r/12)  <  !/ <  24. 


However,  several  properties  of  this  method  must  be  considered.  The  grid 
produced  using  a  simple  shocklike  function  shown  in  Figure  1  is  satisfactory. 
However,  some  minor  skewing  occurs  even  for  this  problem,  where  the  value 
of  A  direction)  was  3  and  the  value  of  the  adaptation  constant  in  the  rj 
direction  was  imity. 

In  Figure  2,  a  grid  is  presented  for  a  rectangular  domain  with  a  circular 
notch  on  the  boundary.  For  this  example,  the  weight  functions  were  set  equal 
to  one  and  the  grid  is  not  adapted  to  any  prescribed  function.  The  results  of 
using  equidistribution  to  construct  a  grid  with  a  unit  weight  function  show 
that  any  discontinuities  in  boundary  slopes  are  propagated  into  the  interior 
grid.  For  comparison,  the  grid  generated  using  Laplace’s  eq^tion 

V^?  =  0,  (17a) 

V*7)  =  0  (17b) 


is  shown  in  Figure  (3).  This  grid  is  very  smooth,  but  does  not  provide  points 
in  certain  regions  of  the  domain.  Comparison  of  the  computed  results  shows 
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that  a  specific  lack  of  smoothness  exists  in  the  independent  arc  equidistribu- 
tion  scheme.  On  the  other  hand,  grid  points  are  uniformly  distributed  on 
constant-{^,  7j)  curves  when  equidistribution  is  used. 

The  last  example  presented  in  this  section  is  that  of  constructing  an 
adaptive  grid  on  a  domain  where  the  adaptation  function  u  goes  from  1  to  0 
along  a  sinusoidal  curve.  From  Figure  4,  it  is  clear  that  significant  grid 
skewness  occurs  in  regions  where  large  adaptation  in  both  directions  is 
requested.  Based  upon  these  examples,  it  does  not  appear  that  the  use  of 
equidistribution  along  independent  arcs  as  a  multidimensional  adaptive  grid 
gf*ne’-'\tor  is  appropriate  (see  Reference  [5]).  Recently,  Eiseman  [9]  has 
presented  results  showing  better  control  of  skewness  by  using  weight  func¬ 
tions  based  on  curvature  of  the  solution  surface.  However,  the  basic  problems 
with  the  equidistribution  scheme  are  not  in  the  choice  of  weight  function  but 
rather  are  due  to  shortcomings  of  the  system  of  governing  partial  differential 
equations. 


POISSON  GRID  GENERATORS 

One  of  the  most  widely  used  grid  generation  schemes  is  based  on  solving 
Poisson  equations  for  the  grid-point  locations  [10].  Winslow  [II]  originally 
employed  Laplace’s  equation  [cf.  Equation  (17)]  to  provide  the  relationship 
between  physical  and  computational  space.  This  approach  provides  the 
smoothest  grid  possible  when  the  smoothness  measure  is  chosen  to  be  the 
integral 


I  =  J  +{^vf]dxdy  (18) 

The  Euler-Lagrange  equations  which  result  from  minimizing  this  integral  are 
the  Laplace  grid  generation  equations  [7]. 

The  grid  structure  produced  using  (17)  is  very  smooth.  A  major  problem 
with  this  scheme  is  that  the  interior  grid  structure  depends  entirely  on 
boundary  data.  Consequently,  it  is  difficult  to  exercise  control  over  the 
interior  point  locations.  Thompson  et  al.  [8]  suggested  that  the  addition  of 
inhomogeneous  terms  (grid  control  functions)  to  (17)  would  provide  a  means 
of  interior  point  control.  In  two  dimensions,  the  Thompson  scheme  is  written 
as 


(19a) 

(19b) 
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where  P  and  Q  are  selected  to  position  the  grid  points  where  desired.  These 
equations  are  usually  solved  for  the  physical  coordinates  by  inverting  the  role 
of  the  dependent  and  independent  variable.  These  equations  are 

=  -  J-[Px^  + Qx^],  (20a) 

ay^  -  =  -J~[Pyi  +  <?y^] .  (20b) 

where 


<x  =  x;  +  y;,  {21a) 

=  +  (21b) 

y  =  x^  +  yf,  (21c) 

and  the  Jacobian  of  the  transformation  is  given  by 

(22) 


Grids  generated  using  (20),  least  in  theory,  can  be  controlled.  This  is 
accomplished  by  appropriate  selection  of  the  grid  control  functions,  P  and 
Q.  One  of  the  problems  with  this  approach  is  to  select  proper  values  of  P  and 
Q  on  the  interior  to  provide  the  desired  controls. 

An  interesting  attempt  at  establishing  the  necessary  values  of  P  and  Q 
was  presented  by  Thomas  and  Middlecoff  [12].  In  their  approach,  the 
boundary  point  distribution  is  assumed  to  be  given,  and  the  boimdary  values 
of  P  and  Q  are  calculated  from  these  given  data.  The  grid  control  functions 
are  then  evaluated  on  the  interior  of  the  domain  by  using  an  appropriate 
interpolation  procedure.  In  order  to  understand  the  details,  it  is  appropriate 
to  review  their  scheme. 

Let  P  and  Q  be  written  in  the  form 


(23a) 

<?=  [T?x  +  Tljl'/'(iT)). 

(23b) 

I 
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With  this  notation,  the  grid  generation  equations  become 

+  +  +  =0,  (24a) 

a[yj£  +  <^>yi]  -2liy^^  +  y[y^^  +  ^Py^]  =0.  (24b) 

The  vector  form  of  these  equations  is  frequently  used.  If 

r  =  ix  +jt/, 

Equations  (24a)  and  (24b)  may  be  written 

a[rjj  +  <i)r^]  -  2/?rj„  +  y  [r„„  +  =0.  (25) 

This  expression  can  be  written  as  a  second-order  differential  equation  for  arc 
length  along  either  constant-^  or  -tj  curves.  To  obtain  the  appropriate 
expression  along  constant-ij  lines,  eliminate  \p  from  (25)  by  taking  the  cross 
product  with  r^.  This  results  in  the  relationship 


“[!/,( +  ■)>»!)  -  +  ■(>!/!)] 


d 

+  y—  — 

1  yj 

<^v\yj 

(26) 


Let  the  local  inclination  of  the  constant-^  curves  be  denoted  by  v  (Figure  5). 
The  local  slope  is  defined  as 


dx 


=  tan  V  = 

£  —  constant 


(27) 


Along  the  constant-^  curves,  the  local  curvature  is 

dv  dv  \ 

dN/d-q' 

That  is,  the  curvature  for  a  given  line  depends  only  on  tj  and  is  directly 
related  to  dv/dt].  From  (27),  this  function  may  be  written 
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y 


If  it  is  assumed  that  the  curvature  of  the  constant-^  lines  is  zero  as  they 
intersect  the  constant-ij  curves,  then  dv/dT)-0.  In  addition,  if  a  local 
orthogonality  is  assumed,  the  right-hand  side  of  (26)  vanishes.  Using  the 
arc-length  derivative  definitions  given  in  (13),  the  final  result  may  be  written 

<pS|  =  0.  (29) 

Along  the  constant-^  curves 

W„  +  +iV,-0.  (30) 

On  the  Tj  =  constant  boundary  the  grid-point  locations  are  specified.  This 
means  that  arc  lengths  are  given.  If  (29)  is  differenced,  the  value  of  (p  may  be 
determined  at  every  point  on  the  boundary.  In  a  similar  manner,  the  values  of 
ip  on  the  constant-^  boundary  may  be  established  from  (30).  Values  for  <p  and 
\{/  on  the  interior  can  now  be  determined  by  using  an  interpolation  scheme. 

ADAPTIVE  GRIDS 

It  is  interesting  to  note  that  the  reduced  equations  (29)  and  (30)  are 
exactly  equidistribution  expressions  given  in  differential  form  as  (11)  and 
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(12).  By  comparing  these  expressions,  the  weight  functions  are  related  to  the 
grid  control  functions  by  the  equations 


or 


1  dw^ 

at  ’ 

(31a) 

1  9wo 

= - ^ 

W2  9t] 

(31b) 

Wj^  = 

(32a) 

Wo  = 

(32b) 

As  noted  previously,  the  concept  of  equidistribution  as  a  grid  generator  has 
not  met  with  great  success  and  is  not  recommended.  However,  the  idea  is  of 
great  value  in  conceptualizing  grid  schemes. 

Since  the  reduced  equations  on  constant-^  and  -tj  surfaces  result  in 
equidistribution  schemes,  it  is  natural  to  wonder  if  the  equations  on  these 
curves  written  without  assumptions  on  curvature  and  orthogonality  might 
look  like  an  equidistribution  scheme.  Without  assumptions  on  the  behavior  of 
the  two  families  of  curves  at  intersections,  these  equations  may  be  written  as 


<#)  —  (/ij  -  2j'^)cot  6  — 


sin  9 


=  0 


(33) 


and 


>/'+(^-2Mjcot^  + 


Sj  sin  9 


=  0. 


(34) 


As  noted  in  Figure  5,  9  is  the  angle  of  intersection  and  fi  and  v  are  the  slopes 
of  the  constant-Tj  and  curves  respectively.  Equations  (33)  and  (34)  may  be 
viewed  as  nonlinear  equidistribution  equations.  The  first  term  in  each  of  the 
square  brackets  is  the  grid  control  function  and  is  related  to  the  original  P 
and  Q  by  (23).  The  second  terms  are  predominately  controlled  by  the 
intersection  angle  of  the  grid  lines.  If  the  grid  lines  intersect  at  90°,  this  term 
vanishes.  However,  the  cot  5  is  multiplied  by  curvature  terms  which  also 


224 


DALE  A.  ANDERSON 


Fig.  6.  Grid  for  sine-wave  problem:  P  and  Q  related  to  u(x,y)  (Equation  33). 


influence  the  magnitude  of  the  expression.  The  last  term  is  primarily  con¬ 
trolled  by  the  curvature. 

If  the  last  two  terms  of  the  bracketed  expressions  in  (33)  and  (34)  are 
small,  the  equidistribution  law  is  essentially  recovered.  Consequently,  if  (p 
and  4'  are  defined  by  (31),  the  grid  becomes  adaptive  if  <|>  and  \p  are 
sufficiently  large.  In  order  to  test  this  idea,  examples  of  grids  computed  using 
the  Thompson  scheme  with  P  and  Q  defined  in  terms  of  a  weight  function 
are  necessary. 

The  first  example  selected  was  the  same  sine  curve  as  previously  used  in 
the  equidistribution  case.  Again,  the  weight  functions  were  selected  to 
depend  on  the  gradient  information  [cf.  Equation  (16)].  The  resiJts  of  using 
Thompson’s  scheme  for  this  case  are  shown  in  Figure  6.  For  large  values  of 
the  clustering  parameters,  the  grid  scheme  produces  a  mesh  that  is  adaptive. 
The  adaptation  is  exactly  as  desired,  following  the  gradient  line  as  prescribed. 
The  shock-wave  grid  of  Figure  1  was  also  recalculated  using  the  adaptive 
scheme  based  on  the  Poisson  generator,  and  the  results  are  shown  in  Figure 
7.  Again,  the  adaptation  to  the  change  in  the  prescribed  function  is  clear.  As 
a  final  example,  the  solution  to  a  more  complicated  adaptation  problem  is 
shown  in  Figure  8.  In  this  case,  the  prescribed  function  was  set  equal  to  zero 
in  the  lower  right  half  of  the  rectangle,  and  equal  to  one-half  on  the  upper  left 
half  of  the  domain  except  for  a  small  semicircular  region  where  the  function 
has  a  value  of  unity.  The  adaptation  to  the  functional  changes  is  remarkable. 
Notice  that  difference  in  the  magnitude  of  the  gradient  across  the  diagonal  is 
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Fig.  7. 


Grid  for  shocklike  function:  P  and  Q  determined  from  u{x,y). 


Fic.  8.  Adaptive  grid  for  bat-shaped  function. 
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clearly  shown  at  the  boundary  of  the  semicircular  region.  The  circular 
boundary  is  also  well  represented,  considering  the  change  from  a  rectangle  to 
a  circular  grid. 

The  Thompson  scheme  can  be  used  in  an  adaptive  manner  as  evidenced 
by  the  examples  calculated.  In  the  solution  of  a  system  of  partial  differential 
equations,  the  control  functions  P  and  Q  are  written  in  terms  of  some 
measure  of  the  solution,  rather  than  a  prescribed  test  function  as  was  done  for 
the  examples  presented  here. 

The  derivation  of  the  relationship  between  the  grid  control  functions  of 
the  Thompson  grid  scheme  and  the  solution  variables  for  an  adaptive  scheme 
was  presented  for  the  two-dimensional  case.  While  the  three-dimensional 
development  is  not  included  here,  the  same  philosophy  leads  to  a  similar 
identification  of  the  grid  control  functions  with  the  weight  fimctions  of  an 
equidistribution  scheme.  The  lead  terms  of  the  Poisson  grid  scheme  written 
along  constant-computational-coordinate  surfaces  are  consistent  with  those  of 
arc  equidistribution. 


SUMMARY  AND  CONCLUSIONS 

The  idea  of  equidistribution  has  been  reviewed  as  a  technique  for  gener¬ 
ating  grids.  The  difficulty  of  using  independent  arc  equidistribution  as  a  grid 
generator  and  selecting  appropriate  weight  functions  was  demonstrated 
through  the  use  of  simple  examples.  Next,  Poisson  grid  generators  were 
reviewed,  emphasizing  the  difficulty  of  selecting  grid  control  functions.  The 
Thompson  grid  generation  equations  were  written  along  constant-computa¬ 
tional-coordinate  surfaces  and  shown  to  be  equivalent  to  a  set  of  nonlinear 
equidistribution  equations.  From  this  point  of  view,  appropriate  choices  for 
the  grid  control  functions  on  the  interior  were  identified,  making  the  stan¬ 
dard  grid  scheme  completely  adaptive.  This  behavior  was  demonstrated  with 
simple  numerical  examples. 

The  concept  of  using  the  Poisson  grid  generators  adaptively  is  significant. 
The  necessary  smoothness  properties  are  part  of  the  grid  generator  (not  the 
weight  function),  and  simple  e.xamples  illustrate  the  success  of  the  approach. 
From  an  economic  view,  most  organizations  do  use  Poisson-type  equations  as 
grid  generators,  and  these  can  be  made  into  solution-adaptive  schemes  with 
near-trivial  changes.  Consequendy,  existing  software  can  be  retained  and  no 
major  investment  in  new  programs  is  necessary  to  obtain  the  ability  to 
produce  adaptive  grids. 
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